ENTROPY  STABLE  APPROXIMATIONS  OF  NAVIER-STOKES  EQUATIONS 
WITH  NO  ARTIFICIAL  NUMERICAL  VISCOSITY 

EITAN  TADMOR  AND  WEIGANG  ZHONG 


Abstract.  We  construct  a  new  family  of  entropy  stable  difference  schemes  which  retain  the  precise 
entropy  decay  of  the  Navier-Stokes  equations,  d/dt  Jx(—pS)dx  =  —  f  ((A +  2 p)qx/6  +  K,(dx/9)2)dx.  To 
this  end  we  employ  the  entropy  conservative  differences  of  [Tadmor2004]  to  discretize  Euler  convective 
fluxes,  and  centered  differences  to  discretize  the  dissipative  fluxes  of  viscosity  and  heat  conduction. 
The  resulting  difference  schemes  contain  no  artificial  numerical  viscosity  in  the  sense  that  their  entropy 
dissipation  is  dictated  solely  by  viscous  and  heat  fluxes.  Numerical  experiments  provide  a  remarkable 
evidence  for  the  different  roles  of  viscosity  and  heat  conduction  in  forming  sharp  monotone  profiles  in 
the  immediate  neighborhoods  of  shocks  and  contacts. 
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1.  Introduction 


We  consider  the  full  Navier-Stokes  (NS) 
mension, 


equations  for  compressible  viscous  flows  in  one-space  di- 
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Here,  p  =  p(x,t )  is  the  density  of  the  flow,  m  =  m(x,t)  is  the  momentum  and  E  =  E(x,t)  stands  for 
the  total  energy  per  unit  volume.  The  three  equations  express,  respectively,  conservation  laws  of  mass, 
momentum  and  total  energy  for  the  flow,  driven  by  convective  fluxes  on  the  left  together  with  viscous 
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and  heat  fluxes  on  the  right.  These  fluxes  involve  the  velocity  q  :=  m/p,  the  pressure  p  =  p(x,  t )  which 
is  determined  by  an  ideal  polytropic  equation  of  state, 

9 

7YI 

P  =  (  7  — l)e,  e:=E-—,  (1.2) 

and  the  absolute  temperature,  0  =  6(x,  t )  >  0,  such  that  Cvp9  =  e.  The  constant  7  >  1  is  the 
specific  heat  ratio  and  e  =  e(x ,  t )  is  the  internal  energy.  On  the  RHS  of  (1.1)  we  have  the  viscous  and 
heat  fluxes,  depending  on  the  constant  Lame  coefficients  of  the  viscosity  A,  p  >  0  and  the  constant 
conductivity  k  >  0.  Finally,  Cv  >  0  is  the  specific  heat  at  constant  volume;  for  simplicity,  we  set 
Cv  =  1  while  rescaling  k  n/Cv. 

If  the  heat  flux  is  excluded  from  the  full  NS  equations,  i.e.  k  =  0,  we  obtain  the  viscous  NS  equations 
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On  the  other  hand,  if  we  turn  off  the  viscosity,  i.e.  A  =  p  =  0,  then  the  NS  equations  amount  to 
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If  the  heat  flux  and  viscosity  are  both  taken  away,  the  system  (1.1)  is  reduced  to  the  compressible 
Euler  equations, 
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The  additional  viscous  and  heat  flux  terms  on  the  RHS  of  the  various  NS  equations  (1.1),  (1.3) 
or  (1.4),  are  dissipative  terms  in  the  sense  that  they  are  responsible  for  the  dissipation  of  the  total 
entropy.  To  this  end,  we  now  discuss  the  entropy  balance  associated  with  the  above  equations.  We 
begin  with  the  specific  entropy  S  :=  In  (pp-7).  A  straightforward  manipulation  on  (1.1),  (1.2)  yields 
the  transport  equation, 

St  +  qSx  = —(In  9)xx  +  (\  +  2p)— +  —  ( .  (1-6) 

p  e  p  V  0  ) 

Multiplied  by  p,  (1.6)  becomes 

pSt  +  mSx  =  n(ln9)xx  + (\  +  2p)^- +  k  (^-^j  .  (1.7) 


On  the  other  hand,  pre-multiplying  the  continuity  equation,  pt  +  mx  =  0,  by  S  and  adding  it  to  (1.7), 

pS^  +  ^(~mS  + ^ln9^x)  =  ~(kX  +  2^<J  ~  '  (1'8) 

Spatial  integration  of  (1.8)  then  yields 

I  (~pS)dx=  -(X  +  ‘2V)  J^l-^dx  -  k  J  9X2  (i)  dx.  (1.9) 

Since  the  expression  on  the  right  is  negative,  we  conclude  that  the  total  entropy,  Jx(—pS)dx,  is 
decreasing  in  time,  thus  recovering  the  second  law  of  thermodynamics,  e.g.,  [GrootMazurl984].  In 
fact,  equation  (1.9)  specifies  the  precise  entropy  decay  rate,  which  is  dictated  by  the  viscous  and  heat 
fluxes  through  their  dependence  on  the  nonnegative  k,  A,  and  p. 


Entropy  stable  approximations  for  the  Navier-Stokes  equations 
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The  purpose  of  this  paper  is  to  present  a  systematic  study  of  difference  schemes  which  respect  the 
above  entropy  dissipation  statements.  The  typical  approach  by  practitioners  in  the  field  of  Computa¬ 
tional  Fluid  Dynamics  is  to  address  the  general  issue  of  entropy  stability  by  adding  ‘enough’  artificial 
numerical  viscosity  —  often  an  excessive  amount  of  it,  in  order  to  mask  various  discretizations  errors 
and  enforce  the  decay  of  the  total  entropy  f  (—pS)  dx.  Our  aim  here  is  to  construct  more  ’faithful’ 
approximations  of  the  NS  equations,  with  a  discrete  analogue  for  the  precise  entropy  decay  statement 
in  (1.9).  Our  prototype  result  is  the  following. 


Theorem  1.1.  Consider  the  NS  equations  (1.1), 

d  d  d 2 

_u+_fH  =  £_d(u), 

Here,  u  =  [p,  m,  i?]T  is  the  vector  of  conservative  variables,  f(u)  is  the  corresponding  3-vector  of 
fluxes  f(u)  =  [m,  qm  +  p,  q(E  +  p)] T ,  and  ed(u)  stands  for  the  combined  viscous  and  heat  fluxes, 

ed(u)  =  (A  +  2 p)  [0,  q,  q2 / 2] T  +  k  [0, 0,  9 ] T 

where  e  signals  the  vanishing  amplitudes  of  viscosity  and  heat  conductivity.  We  approximate  these  NS 
equations  by  a  semi- discrete  scheme  of  the  form 
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Let  f*  i  =  f*(uj,,  u„+i)  is  the  numerical  flux  given  the  by  the  explicit  formula, 
2 


f;+4 = u  - 1)  e 


mP+1  — 


— (  (£J,v„+1  -  v„) 


-fj. 


v  =  v  u  :  = 


-|-s+7  +  i,|,-1]t.  (i.iob) 


Here,  [£J  =  ^+1}|=1  are  three  linearly  independent  directions  in  v -space  at  our  disposal  (consult 

examples  3.1,  3.2  and  3.3  below);  {r7  =  r^+1  }|=1  is  the  corresponding  orthogonal  system  and  {m^  = 

mi+±}j=i  are  intermediate  values  of  the  momentum  specified  along  the  corresponding  path,  vJI+1  = 

v-1  +  (F,vy+i  —  starting  with  v1  =v^  and  ending  with  v4  =  vi,+i.  Then,  the  resulting  scheme 

(1.10a),  (1.10b)  is  entropy  stable  and  the  following  discrete  entropy  balance  holds 1 


d_ 

dt 


'^j(-pvSu)Axv  = 


-(A  +  2/j)  ^ 

V 


'A%+± 


Ax"H 


(l/»)(,4vj<°.  (l'H) 


We  briefly  describe  the  content  of  the  paper,  leading  to  the  above  result.  The  entropy  balance  (1.11) 
is  a  precise  discrete  analogue  of  (1.9).  The  scheme  (1.10a), (1.10b)  contains  no  artificial  numerical 
viscosity  in  the  sense  that  entropy  dissipation  is  driven  solely  by  the  viscous  and  heat  fluxes.  In  the 
particular  case  that  viscosity  and  the  heat  conduction  are  absent,  n  =  \  =  p  =  0,  then  the  entropy 
balance  (1.8)  is  reduced  to  the  formal  entropy  equality  of  the  Euler  equations  (1.5), 

§iCps)  +  m(-mS)=0’  <L12> 


1We  let  z,,,  i  and  z„,  i  denote  the  arithmetic  and  geometric  means,  i 
2  '  2  '  2 


{zv  +  z„+ 1)/2  and  zv+ 1  =  yJzvzv+\. 
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which  in  turn,  implies  the  entropy  conservation  Jx(-pS)(-,t)  dx  =  f  (— pS)(-,  0)  dx.  Similarly,  setting 
ed  =  0  we  omit  the  dissipative  terms  in  NS  equations,  and  the  difference  scheme  (1.10a)  becomes 
entropy  conservative , 

^2{~PuSu){t)  Axv  =  y^X~puS„)(0)Axu. 

Entropy  conservative  schemes  are  studied  in  section  3,  following  [Tadmor2004] .  The  key  ingredient  here 
is  the  construction  of  their  entropy  conservative  fluxes,  such  as  f*+1  in  (1.10b).  These  fluxes  employ 

the  so  called  entropy  variables,  v  =  v(u),  which  are  discussed  in  section  2.  The  main  results  are  then 
summarized  in  theorems  3.1  and  3.2.  Finally,  in  section  4  we  present  a  series  of  numerical  simulations 
with  the  new  schemes.  The  entropy  conservative  approximations  of  Euler  equations  are  ‘purely  dis¬ 
persive’  and  as  such,  their  solutions  experience  dispersive  oscillations,  interesting  for  their  own  sake, 
consult  [Laxl986,  HFMM1986,  Tadmorl986,  LLV1993,  LevermoreLiul996,  LeFlochRhode2000]  and 
the  references  therein.  Turning  to  the  NS  equations,  our  simulations  provide  a  remarkable  evidence 
for  the  different  roles  that  viscosity  and  heat  conduction  have  in  removing  the  dispersive  oscillations, 
to  yield  sharp  monotone  profiles  of  well-resolved  shock  and  contact  layers.  No  limiters  were  added, 
but  instead,  the  viscous  and  heat  conduction  terms  in  NS  equations  are  found  to  serve  as  accurate 
edge  detectors 

Remark  1.1.  The  viscous  NS  equations  dissipate  a  general  family  of  entropies,  —ph(S),  where  h(S) 
is  an  arbitrary  increasing  function.  Indeed,  arguing  along  the  above  lines  we  multiply  the  continuity 
equation  by  h(S)  and  adding  it  to  (1.7)  x  h'(S)  to  find 

r\  r\ 

—  (-  ph(S))  +  —^-mh(S)  +  K(lnd)xh'(S)^  = 

=  Kh"(S)Sx^-h’(S )  ^(A  +  2m)|  +  n(^^  .  (1.13) 

In  the  case  that  the  heat  conduction  is  absent,  the  first  term  on  the  RHS  of  (1.13)  vanishes,  and  we 
are  left  with 

Jtf  i~  ph(S))dx  =  -(A  +  2  p)  J  | h\S)dx .  (1.14) 

Thus,  the  viscous  NS  equations  imply  the  dissipation  of  a  family  of  entropies,  J  (-ph(S))  dx  for  all 
h\S)  >  0;  consult  [Hartenl983].  Each  one  of  these  entropies  carries  its  own  entropy  conservative 
flux  f  *  i .  The  explicit  construction  of  such  fluxes  is  outlined  in  theorem  3.1  below.  Combining 

these  entropy  conservative  fluxes  together  with  centered  differencing  of  the  additional  viscous  terms, 
(A  +  2/x)[0,  q ,  q2 / 2]T,  yield  a  generalization  of  theorem  1.1  which  recovers  the  precise  entropy  balance 

(  hi  1). 

We  note  in  passing  that  when  heat  conduction  is  present,  however,  the  negativity  of  the  first  term 
on  the  right  of  (1.13)  requires  h”(S)  =  0,  so  that  we  are  left  with  one  canonical  entropy,  h(S)  ~  S 
discussed  in  theorem  1.1;  consult  [HFMM1986]. 

Remark  1.2.  It  is  straightforward  to  generalize  the  recipe  for  ‘faithful’  entropy  stable  approximations 
of  multidimensional  NS  equations.  The  extension  is  carried  out  dimension  by  dimension  and  as 
indicated  in  the  one-dimensional  setup  of  theorem  1.1,  one  has  the  freedom  of  choosing  different  paths 
in  phase  space. 


2.  Entropy  dissipation 


2.1.  The  entropy  variables.  We  turn  our  attention  to  general  systems  of  hyperbolic  conservation 
laws 


d  d 

-u+-f(u)  =  0. 


(x,  t)  6  R  x  [0,  oo), 


(2.1) 
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governing  the  IV-vector  of  conserved  variables,  u  =  [u\, . . .  ,  uat]t.  The  entropy  equality 

!"<“)  +  lF(u)  =  °-  <2-2> 

is  an  additional  conservation  law  for  an  entropy  function,  U( u),  which  is  balanced  by  an  entropy  flux 
F( u).  The  Euler  equations  (1.5)  are  viewed  as  a  prototype  example  for  such  systems,  with  the  three 
conservative  variables  u  =  [p,  to,  E] T  balanced  by  the  flux  f  =  [to,  qrn  +  p.  q(E  +  p)] T  and  endowed 
with  entropy  pairs  {U,  F)  =  (  —  ph(S),  — mh(S )).  We  now  briefly  recall  the  circle  of  ideas  linking  the 
dissipation  of  the  total  entropy,  f  U(u(-,t))dx,  and  the  realization  of  u  as  a  vanishing  viscosity  limit, 
in  analogy  to  the  vanishing  NS  limits  and  their  relation  to  entropic  solutions  of  the  Euler  equations. 
We  refer  to  e.g.,  [Liul991]  and  [Dafermos2000],  for  a  more  comprehensive  discussion. 

Let  (U(u),  F( u))  be  a  given  entropy  pair  associated  with  (2.1).  Note  that  U(u)  satisfies  the  entropy 
equality  (2.2)  if  and  only  if  it  is  linked  to  an  entropy  flux  function  E( u)  through  the  compatibility 
relation 

^fu  =  EuT.  (2.3) 

Indeed,  multiplying  (2.1)  by  on  the  left,  one  recovers  the  equivalence  between  (2.2)  and  (2.3)  for 
all  u’s  solving  (2.1).  Of  course,  these  formal  manipulations  are  valid  only  under  the  smooth  regime. 
To  justify  these  steps  in  the  presence  of  shock  discontinuities,  the  conservation  law  (2.1)  is  realized 
by  appropriate  vanishing  viscosity  limits.  To  this  end,  we  define  the  entropy  variables  v(u)  :=  Hu( u). 
We  make  the  additional  assumption  that  the  entropy  U  (u)  is  convex ,  so  that  the  mapping  ue>v  is 
a  one-to-one.  Following,  [Godunovl961,  Mockl980],  we  claim  that  the  change  of  variables,  u  =  u(v), 
puts  the  system  (2.1)  into  the  equivalent  symmetric  form, 

iu(v)  +  l;f(u(v))  =  °-  <2-4) 

The  system  (2.4)  is  symmetric  in  the  sense  that  the  Jacobian  matrices  of  fluxes  are2, 

uv(v)  =  (uv(v))T  and  fv(v)  =  (fv(v))T.  (2.5) 

Indeed,  a  straightforward  computation  utilizing  the  compatibility  relation  (2.3),  shows  that  u(v)  and 
f(v)  are,  respectively,  the  gradients  of  the  corresponding  potential  functions,  (p  and  ip, 

u(v)  =  </>v(v),  <Mv)  :=  (v, u(v))  -  U (u(v)),  (2.6) 

f(v)  =  V’v(v),  V’(v)  :=  (v,f(v))  -  E(u(v)).  (2.7) 

Hence  the  Jacobian  matrices  H(v)  :=  uv(v)  and  B(y)  :=  fv(v)  in  (2.5)  are  symmetric,  being  Hessians 

of  the  potentials  </>(v)  and  ip{v).  Moreover,  the  convexity  of  U(-)  implies  that  H  is  positive  definite, 

H  =  (Huu)-1  >0. 

Physically  relevant  solutions  of  (2.1)  are  postulated  as  limits  of  the  vanishing  viscosity  solutions, 

!u‘ +  sf(u<)  =  e^d(u,)’  <2-8) 

where  d(u)  is  any  admissible  dissipative  flux,  and  e  j  0  stands  for  vanishing  amplitudes  such  as  the 
viscosity  coefficients  A,  p,  the  heat  conductivity  k,  etc.  Here,  the  admissibility  of  the  dissipative  flux 
requires  the  Jacobian  du  to  be  H-symmetric  positive-definite,  that  is, 

d h  =  (d h)T  >  0,  dH  :=  d UH.  (2.9) 

2For  brevity  of  notations  we  often  write  f(v)  for  f(u(v))  whenever  the  different  dependence  of  f(u)  and  f(v)  is  made 
clear  by  the  distinction  between  the  conservative  variables  u  and  entropy  variables  v. 
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If  we  express  the  dissipation  flux  in  terms  of  the  entropy  variables,  d(v)  =  d(u(v)),  then  admissi¬ 
bility  requires  that  the  v-Jacobian  of  this  flux  will  be  positive  symmetric,  d//  =  dv(v)  =  d^(v)  >  0. 
Thus,  in  this  case  (2.8)  reads 


lu(v<) +  Kf(v,)  =  e^d(v!)' 

and  all  v-dependent  fluxes  —  the  temporal,  spatial  and  the  dissipation  flux  have  symmetric  Jacobains. 

We  now  integrate  (2.8)  against  vT  =  C/J”,  employ  the  compatibility  relation  fx  =  Fj  and  use 
‘differentiation  by  parts’  on  the  admissible  dissipation  on  the  right  to  find 

^(ue)  +  ^(^K)  “  e(v£’d(ue)*>)  =  -<Kvz,  dvv*)  <  0.  (2.10) 

Letting  e  j  0,  we  obtain  the  entropy  inequality,  [Laxl973] 

!;/(“) +  L.F(u)<„,  (2.H) 

which  is  the  generalization  of  the  entropy  decay  statements  for  NS  equations  (1.8)  and  (1.13). 


2.2.  Examples  of  entropy  pairs  for  Navier-Stokes  equations.  The  NS  equations  admit  the 
family  of  convex  entropy  pairs 

U (u)  =  -ph(S),  F{ u)  =  -mh(S),  h'(-)  >  0.  (2.12) 

Here  S  =  In  (pp-7)  is  the  specific  entropy  and  the  convexity  of  the  corresponding  U( u)’s  as  functions 
of  u  =  ( p,m,E)T  holds  iff  h'(S )  —  7 h"(S)  >  0,  [Hartenl983].  We  consider  two  prototype  examples. 

Example  2.1.  The  simplest  choice  of  h(S)  is  the  specific  entropy  S  itself, 

h(S)  =  S  =  In  (pp"7)  .  (2.13) 

Straightforward  computation  gives  us  the  following  entropy  pair,  entropy  variables,  and  potentials. 

•  Entropy  pair  U( u)  =  —  pS  and  F( u)  =  — mS ; 

•  Entropy  variable  (consult  [Hartenl983]) 


v(u)  = 


-E/e  -5  +  7  +  1 
q/e 
-i/e 


(2.14) 


with  the  inverse  mapping 

p 


u  v  = 


7-1 


-v3 

—  V3 

V2 

=  w 

V2 

1  v2 

L  2v3  J 

1  v2 

L  1  2v3  J 

,  where  w  = 


_  ( 7-1  y-1^ 


(-^3)7 


c  ,  v2 

,  S  =  7  -  Vi  +  — ; 

2v3 


•  Potential  pair  (j>  =  ( 7  —  l)p  and  ip  =  ( 7  —  l)m. 

In  this  case,  the  general  statement  of  entropy  balance  in  (2.10)  with  the  entropy  pair  (U,  F)  = 
(— pS ,  — mS )  amounts  to  the  one  we  have  in  (1.9), 


j  -{pS)dx  =  -e  J  (vex,  d Vvex)dx  =  -(A +  2 p)  J  j- 


2  r  a2 

dx  —  k  /  dx  <  0.  (2-15) 


Example  2.2.  A  particularly  convenient  form  of  entropy  variables  is  associated  the  entropy  function 
(consult  [Hartenl983,  Tadmor2004] ) , 


U(u) 


— ph(S )  with  h(S) 


7  +  1 
7-1 


s 

e7+i  ? 


(2.16) 


where  we  have  the  following. 
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1  '"j  l  ]_  ~\~  /"y  l 

Entropy  pair  U( u)  =  - - (pp)  x+r  and  F(u)  =  - - q(pp)1+j', 

—  7 

Entropy  variable  v(u)  =  Vut/(u)  =  —(pp)  1+r 
with  the  inverse  mapping 

7 

u(v)  =  —  (pp)r+ 1 


1-7' 

E 

—m 

P 


^3 

-Vi 

—  _ 

(7  - 1)  1 

1 

1 

CO 

7 

1-7 

vs 

-V2 

Vi 

\  /  j 

u  Vl 

•  Potential  pair  ((f),  if))  =  ((pp)  T+1 ,  m(pp  7)r+i 
In  case  that  the  heat  conduction  is  absent  (k  =  0),  we  apply  the  general  statement  of  entropy  balance 
(2.10)  with  the  entropy  pair,  (U,  F)  =  ]r^((pp)~ ,  q(pp)JF/),  obtaining 


d_ 

dt 


~(pp) 1+7  dx  =  -e  /  (vex,  d Vvex)dx 

J  X 

(7-l  )h'(S)ql 


=  “ +  2 p)  f 

J  X 


(1  +  7  )0 


dx  =  — 


\  +  2p 

1  +  7  , 


eS/a+TO^ 


dx  <  0. 


(2.17) 


Remark  2.1.  As  noted  in  [Hartenl983],  the  flux  f(v)  is  a  homogeneous  function  of  degree  r/  =: 
(1  +  7)/(l  —7),  f(av)  =  a'?f(v),  Va  G  M.  Homogeneity  implies  that  fv(v)v  =  r] f(v)  which  in  turn, 
enables  us  to  rewrite  the  spatial  flux  in  (2.1)  in  a  skew-adjoint  form,  f(u)j;  =  (fvvx  +  (fvv)x)/(r?  +  l); 
consult  [Tadmorl984a] . 


3.  Entropy  stable  schemes  for  Navier-Stokes  equations 


3.1.  Entropy  conservative  schemes.  We  turn  our  attention  to  consistent  approximations  of  (2.1), 
based  on  senri-discrete  conservative  schemes  of  the  form 


d  ,  . 

dtuJt) 


(3.1) 


Here,  u„(f)  denotes  the  discrete  solution  along  the  grid  line  (xv,t),  Axu  :=  \(xv+\  —  xu-\)  is  the 
possibly  variable  mesh  spacing  and  fj,+  i  is  the  Lipschitz-continuous  numerical  flux  which  occupies  a 
stencil  of  (2 p  +  l)-gridvalues, 


—  f  (W— p+l)  J  W+p)- 

The  scheme  is  said  to  be  consistent  with  the  system  (2.1)  if  f  satisfies  f(u,  u,  •  •  •  ,  u)  =  f(u),  Vu  G  Rn. 
By  making  the  changes  of  variables  =  u(vi/),  we  obtain  the  equivalent  form  of  (3.1) 

Ttu{Mt)]  =  -lk(Fi-f-i)'  <3-2) 

The  essential  difference  lies  with  the  numerical  flux,  f;y+ 1 ,  which  is  now  expressed  in  terms  of  the 
entropy  variables, 


f^+i  =  f  K-p+i,  •  •  •  ,  v^+p)  :=  f  (u  (v^_p+i) ,  •  •  •  ,  u  (v^+p))  , 
consistent  with  the  differential  flux, 

f(v,v,...  ,v)  =  f(v)  =  f(u(v)). 


(3.3) 


The  senri-discrete  schemes  (3.1)  and  (3.2)  are  completely  identical.  It  proved  useful,  however,  to  work 
with  the  entropy  variables  rather  than  the  usual  conservative  ones,  since  system  (2.1)  is  symmetrized 
with  respect  to  these  entropy  variables.  The  entropy  variables-based  formula  (3.2)  has  the  advantage 
that  it  provides  a  natural  ordering  of  symmetric  matrices,  which  in  turn  enables  us  to  compare  the 
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numerical  viscosities  of  different  schemes,  consult  [Tadmorl984b,  Tadmorl987]  for  details.  In  partic¬ 
ular,  we  will  be  able  to  utilize  the  entropy  conservative  discretization  of  [Tadmor2004]  for  the  Euler 
convective  part  of  the  equations,  and  thus  recover  the  precise  entropy  balance  dictated  by  viscous  and 
heat  fluxes  in  NS  equations. 

Let  (U,  F)  be  a  given  entropy  pair.  We  proceed  with  the  construction  of  an  entropy  conservative 
scheme,  in  the  sense  of  satisfying  a  discrete  entropy  equality  analogous  to  (2.2), 

i.U(Mt))+A-^+i-K_i)=0.  (3,4) 

Here,  Fu+i  =  F  (u„_p_|_i,  •  •  •  ,  u^+p)  is  a  consistent  numerical  entropy  flux,  such  that  F( u,  u,  •  •  •  ,  u)  = 

.F(u),  Vu  £  Rw.  The  numerical  flux  of  such  entropy  conservative  schemes  will  play  an  essential  role 
in  the  construction  of  entropy  stable  schemes,  by  adding  a  judicious  amount  of  physical  viscosity. 

The  key  step  in  the  construction  of  entropy  conservative  schemes  is  the  choice  of  an  arbitrary 
piecewise-constant  path  in  phase  space,  connecting  two  neighboring  gridvalues  v„  and  v„+i  through 
the  intermediate  states  {v^+i  at  the  spatial  cell  [x„,  xv+i).  Let  }|Li  be  an  arbitrary  set  of  N 

linearly  independent  JY-vectors,  and  let  {i3  ,  }f=1  be  the  corresponding  orthogonal  set.  We  introduce 

the  intermediate  gridvalues  {v^+1}j^.1,  which  define  a  piecewise  constant  path  in  phase  space 


v1  , 
U+2 


=  v,„ 


7  +  1  j 

v  ,  !  =  V 


v+i 
N+ 1 


.  +  ( e. 


=  v„+i 


Yj  4  j  =  1, 2,  • 


,1V  —  1,  Av 


^+i  :=  v^+i  -  v„. 


(3.5) 


Theorem  3.1  (Tadmor2004,  Theorem  6.1).  Consider  the  system  of  conservation  laws  (2.1).  Given 
the  entropy  pair  (U,  F),  then  the  conservative  scheme 


with  a  numerical  flux  f* ,  i 


d  ,  . 


v+ 


1 

2 


N  ^ 

E- 


3  = 1 


(3.6) 


(3.7) 


is  an  entropy -conservative  approximation,  consistent  with  (2.1), (2. 2).  Here,  v  are  the  entropy  vari¬ 
ables,  v  =  f7u(u)  and  ip(v)  is  the  entropy  potential  (2.7)  ip(v)  =  (v,f(u(v)))  —  F(u(v)). 


The  proof  is  based  on  the  requirement  of  entropy  conservation  in  [Tadmorl987,  Theorem  5.2], 
(Av„+1  ,r+i)  =  A^+1,  A^+1  :=ip(vv+1)  -  ^(v„). 

The  numerical  flux  (3.7)  satisfies  this  entropy  conservation  requirement,  for 


Av^’C+i 


n 

E- 

3= 1 
N 

3= 1 


vS )  -  ^ 


VJ  ! 

"+3 


i  >  Av 


u+ 


Y+1 


v+i 


f3  ,  i ,  Av,  ,  i 


I/+7 


-Ip 


=  1l> 


- 

"+1 J 


v'W  I  -  ip  ( v ' 


i/+7 


=  AVw±- 


(3.8) 
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In  addition,  f*  ,  is  a  consistent  flux  satisfying  (3.3).  Indeed,  if  we  let  v'  \ (£)  denote  intermediate 


v+  ^ 

path,  v£?tf)  :=  (vJ 


12+ i 


+  v^11)/2  +  e<Ci,Av_i) 


u+i 


12+ 4 


we  have 


-§^r 


connecting  and  v^1, ,  then  by  (2.7), 


1+1 
v  1 
U+2 


VJ  ! 
V+2 


J KO  # 


U£>  #,r 


Inserted  into  (3.7),  we  can  rewrite  the  entropy-conservative  flux  (3.7)  in  the  equivalent  form 


f*  i  = 


V+ 


t=-h 


'J  +  |(0  )  d£,r3 


i3 


(3.9) 


N 


and  the  consistency  relation  (3.3)  now  follows,  f*(v,  v)  =  (  f(v) 


i=i 


=  f(v)-  D 


We  emphasize  that  the  recipe  for  construction  entropy-conservative  fluxes  in  (3.7)  allows  an  arbitrary 
choice  of  a  path  in  phase  space.  We  demonstrate  this  recipe  with  three  examples. 


Example  3.1.  Set  {rJ}  along  the  standard  Cartesian  coordinates,  r7  1  =  ej,  j  =  1, 2, . . .  ,  N.  In 


this  case  we  have 


VJ  x  = 
U+2 


Mi>-  ,K) 


lT 


N 


j  =  2,3,...  ,1V-  1, 


and  the  entropy  conservative  flux  (3.7)  is  given  by  the  particularly  simple  explicit  formula 


f*  i  = 


^(vLi)  -  VM) 


12+ i 


) 


V>(v„+1)  -^(v^+i) 


-1  T 


K+i)1-(v,)1  ’  (v,+1) ,-K). 


(v^+Ojv  -  (Vl/) 


N 


(3.10) 


We  carried  out  numerical  experiments  with  these  fluxes  for  the  approximate  solution  of  Euler  and  NS 
equations.  The  formulation  is  particularly  simple  though  the  computed  intermediate  values  might  lie 
outside  the  physical  space  p,p  >  0. 

Example  3.2.  A  more  ‘physically  relevant’  choice  than  the  Cartesian  path  is  offered  by  a  Riemann 
path  which  consists  of  {iE  +  l  }fE1,  stationed  along  an  (approximate)  set  of  right  eigenvectors,  {r^+1 }, 

of  the  Jacobian  fu(u„,i ).  Set  vJ  ,  =  v(iE  7  =  1,2, 

'  u+2  12+  ^  V  u+^ 

This  will  be  our  choice  of  a  path  for  computing  entropy  stable  approximations 


,  N,  and  let  ’s  be  the  orthogonal  system 

"  1  2  •'To  ‘'To 

to  {v^+1  - 

of  NS  equations  in  section  3.2  below.  The  resulting  flux,  mixing  conservative  and  entropy  variables 
admits  the  somewhat  simpler  form 


N  E 

f;+i  =  E- 

3= 1 


lE+\ 

12+b 


-  E 


lE  ] 

12+k 


£3  i  where  E(u)  =  (l/u(u))Tf(u)  —  F( u) 


l3  i ,  Av 


12+ 


12+i 


12+ 


(3.11) 


Example  3.3.  If  all  rJ’s  are  chosen  to  approach  the  same  direction  of  Av^i,  then  by  (3.9)  the  flux 
(3.9)  ‘collapses’  to  the  entropy-conservative  flux 


f,  i  = 


'£=- 


Vri(£))  d£, 


v„+i(0  :=  g (v„  +  viH_i)  +  £Avi/+i. 


(3.12) 
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The  resulting  flux  (3.12)  was  introduced  in  [Tadmorl986]  and  was  the  forerunner  for  the  family  of 
entropy  conservative  fluxes  outlined  in  theorem  3.1.  It  has  the  drawback,  however,  that  its  evaluation 
requires  a  nonlinear  integration  in  phase  space.  Thus,  with  the  loss  of  linear  independence,  we  lose 
here  the  explicit  evaluation  of  the  entropy  conservative  flux  offered  in  (3.7)  and  demonstrated  in  the 
previous  two  examples. 

3.2.  Entropy  stable  semi-discrete  schemes  for  Navier-Stokes  equations. 

3.2.1.  The  compressible  Euler  equations.  Let  (U,  F )  be  an  admissible  entropy  pair  associated  with  the 
Euler  equations  (1.5),  let  v  =  v(u)  denote  the  corresponding  entropy  variables  outlined  in  examples 
2.1  and  2.2  above.  To  conserve  the  total  entropy  f  U(u(-,t))dx,  we  appeal  to  the  semi-discrete  scheme 
(3.6)  with  the  entropy-conservative  numerical  flux  (3.7), 


3 


E 


To  compute  f *  i ,  we  distinguish  between  two  cases.  If  wv  =  v„+i,  we  employ  the  equivalent  form  of 
the  numerical  flux  in  (3.9), 


v+ 


1 

2 


3 


E 


which  implies  that  all  the  intermediate  gridvalues  coincide,  v„  =  v1  ,  =  v2  ,  =  v3  ,  =  v 1  ,  =  v,,+i 

V+±  V+±  V+±  T 

and  the  entropy-conservative  flux  amounts  to  f*  i  =  f„  =  f,y+ 1 .  Otherwise,  if  \v  /  v^i,  we  choose 

to  work  along  the  path  which  is  dictated  by  an  (approximate)  Riemann  solver.  Specifically,  we  use 
the  eigensystem  of  the  Roe  matrix,  [A^i  =  -A(u„,  u„+i),  [Roel981], 


eE-o3 
2  Vi 


0 

In In2 
2  %+l 


1 


(3  -  t)%+\ 
Hu+i  Hu+i  -  (7  -  1  )o2 


0 

7-1 


j 


(3.13a) 


Here  q  and  H  are  the  average  values  of  the  velocity  q  and  total  enthalpy  H  =  ( E+p)/p  at  Roe-average 
state, 


%+h  = 


Quy/Pu  +  Qv+ly/  Pu+ 1 
\[Pv  +  y/Pv+ 1 


H^  = 


Huy/ Pv  T  Hu+  \y/ Pu+1 

yfPv  +  y/Pu+l 


(3.13b) 


The  rJ’s  are  the  right  eigenvectors  {r-7  =  rJ  i}?=i  of  the  Roe  matrix  (3.13a)  given  by  (omitting  the 

i2+ 2  2 


subscript  (•)  ,i  of  all  averaged  variables) 


r  1 

r  i  i 

r  1 

r1  = 

q  —  c 

H  —  qc 

,?2  = 

i 

CM 

1^ 

_ i 

,  r3  = 

q+c 

H  +  qc 

with  the  corresponding  left  eigenvector  set  {£J  =  given  by 

2  J 


It  = 

\  (2  +  5)/(4(7-l))  1 

/s.9 

-  l-52/(2(7-l))  ' 

"  -(2-5)/(4(7-l))  1 

-(l  +  5)/{2c) 

X  = 

5/c 

,£  = 

(1  -  8) /(2c) 

(7-1)/(2c2) 

-(7-l)/c2 

(7-1)/(2c2) 

(3.13c) 


(3.13d) 
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Here  5  :=  (7  —  1  )q/c,  and  c  is  the  average  sound  speed  given  by  c2  =  (7  —  1)  ^ H  — 

We  are  now  able  to  form  the  intermediate  path  in  u— space  as  in  (3.5) 

uj+1  =  uj  +  (1\  AuJlP,  j  =  1,2,3.  (3.14) 

Since  the  mapping  between  u  and  v  is  one-to-one,  then  these  intermediate  gridvalues  in  u— space, 
{u?}j=1,  correspond  to  intermediate  gridvalues  {v-?}j=1  in  v— space.  We  let  {rJ}|=1  be  the  (right) 
vectors  connecting  these  v- values,  r-1  :=  vJ+1  —  v-7,  and  let  {£3}j=1  be  the  corresponding  (left)  or¬ 
thogonal  set.  We  summarize  the  algorithm  of  computing  the  entropy-conservative  flux  f*  1  in  the 

v+2 

following. 

Algorithm  3.1.  If  u„  =  u„+i  then  f*+1  =  f(v^);  else 

•  Set  u*+1  :=  \iu  and  compute  recursively  the  intermediate  states, 


1  >  r 


; 3 


Ui+\  =  UJ  x  ,  1  i  Au  ,  i  /  *  1  • 

V+\  V+±  \  "+2’  u+2  /  V+\‘ 


1  •  j  —  1)2,3 

2 

Here,  {^+i }  and  {r^+1  }  are  the  left  and  right  eigensystems  of  the 
Set  r'  j  =  v(uJ+11 )  —  v( 

i/4-  —  v  1/4-  —  '  v 


(3.15) 

Roe  matrix  in  (3.13c),  (3.13d). 


=  J 


’jki 


W  +  1)  and  compute  {£J}^=1  as  the  corresponding  orthogonal  system, 

(3.16) 

"  1  2 

,  (  1+1  \  ,fi\ 

3 

Compute  the  entropy-conservative  numerical  flux,  f*  1  = 

3= 1 


ri  •  yi  •  1  _  vi 


3  “  V’ 


ej 


Av 


^+0 


Remark  3.1.  Observe  that  if  {u^1  —  U'J}^=1  are  linearly  independent  then,  since  uv  is  symmetric 
positive  definite,  the  corresponding  set  of  directions  in  v-phase  space,  {v7+1  —  vJ  }j=  j ,  are  also  linearly 
independent,  at  least  when  uv+i  is  in  a  small  neighborhood  of  u„.  It  guarantees  the  existence  of  the 
orthogonal  set  But  what  happens  when  (^,Au)  =  0  for  certain  f  s?  for  example,  if  u„  is 

connected  to  u^+i  through  a  fc-shock,  then  the  Roe  matrix  [A]^+i  retains  the  perfect  resolution  of 

such  a  shock  by  enforcing  (£  ,  Au)  =  0,  \/j^k  and  we  can  omit  the  contribution  of  these  sub-paths  to 
the  conservative  flux  f*.  The  general  approach  is  to  construct  a  precise  mirror  image  of  the  Roe-path 
in  v-phase  space  in  terms  of  the  right  and  left  orthogonal  systems, 


1 V  :=  [H\u+1_£j ,  j  =  1,2,...,  (3.17) 

2  2 

where  [H] u+i  denotes  an  averaged  symmetrizer  such  that  Au^i  =  [H]v+i  AvJ/+i  (and  there  are 

many  such  averages).  Then,  {r-?}^=1  forms  the  path  in  v-phase  space,  vJ+1  =  v-7  +  (£J ,  Av)rJ",  which 
retains  the  desired  Roe  property  of  perfect  resolution  of  shocks.  Indeed,  if  Au  is  a  fc-shock  with 
speed  s  then  it  satisfies  (omitting  subscripts)  Af  =  sAu  =  s[E]Av.  But  the  Roe  matrix  in  (3.13a) 
is  constructed  so  that  [Roel981],  Af  =  [A]Au  =  [A][H]Av,  and  we  conclude  Av  =  rfc.  Thus, 
(F,Av)  =  0,  Vj  7^  k.  The  corresponding  entropy  conservative  numerical  flux  reads 


{3 


E 

I 


V’ 


+ 


v+i 


e  = 


Av 


v+h 
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3.2.2.  The  Navier-Stokes  equations.  We  turn  to  the  construction  of  entropy-stable  schemes  for  the  full 
NS  equations  (1.1).  To  this  end,  we  rewrite  the  equation  as  a  system  of  conservation  laws 


d  d  d2  9  m 

mu+  oTf(u)  =  eWTd(u)>  u=  rn  >fH  =  (l'm  +  P 
dt  dx  dx  E  q{E  +  p) 


(3.18a) 


with  additional  diffusive  terms 


0  0 

ed(u)  :=  (A  +  2fi)  q  -\- k  0  (3.18b) 

.  «2/2  J  L  6  . 

For  the  convection  part  on  the  LHS,  we  use  the  same  entropy-conservative  differencing  used  for  the 
Euler  equations.  For  the  dissipative  terms  on  the  RHS,  we  employ  standard  centered  differences.  We 
arrive  at  our  main  result. 


Theorem  3.2.  Let  (U,  F )  be  a  given  entropy  pair  of  the  NS  equations  (3. 18a),  (3. 18b),  which  respect 
the  entropy  inequality  (2.10).  Consider  the  semi-discrete  approximation 


jtMt)  +  aL;  (fA  "-*)■  i  ( 


AxV+i  AVi 


Here  f*  i  is  an  entropy  conservative  numerical  flux  (3.7), 


f*  !  = 


3  if  ,  -  l/)[  vf  j 


^+I’AV-+| 


which  is  outlined  in  algorithm  3.1  above. 

{*}  The  resulting  scheme  (3. 19a), (3. 19b)  is  entropy-dissipative  in  the  sense  that 


Y^Uflt)  Ax.  = 


e  /  Ad,+i  \ 

- (  Av  i ,  --^Av  ,i  )  <  0. 

wA  +5  Aw+I  +v 


(3.19a) 


(3.19b) 


(3.20) 


This  entropy  balance  is  a  discrete  analogue  of  the  entropy  balance  statements  (2.15)  and  (2.17). 

{ii}  In  the  specific  case  of  the  canonical  entropy  pair  (U,  F)  =  (— pS ,  —mS),  the  entropy  decay 
(3.20)  amounts  to  (1.11) 

—  Uflt)  Ax.  = 


.1  A*  m 


AX»+I 


,Ax  m<0.  (3.21) 


Proof.  We  multiply  (3.19a)  by  [Uu\ X  =  vj,  then  sum  up  all  spatial  cells  to  get  the  balance  of  the 
total  entropy, 


-  ]T  uv(t)  Ax.  +  ]T  (v„,  r+ 1  -  f 


-  t .  i  )  =  e 


d.+i  -  d.  _  d.  -  d._i 

Ax^i  Axu_i 


(3.22) 


Since  we  chose  f *  x  as  the  entropy  conservative  flux,  a  straightforward  manipulation  on  the  entropy 
^+2 

conservation  requirement  (3.8)  yields  the  conservative  difference, 


V.,f*  x  -f*  ! 

u+i 


=  Fv+i-Fv_i, 


(3.23) 
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where  2Fu+i  =  ^(v„  +Vj,_|_i)  ,  f„+iy  —  +  if)(\ru+ 1)).  On  the  other  hand,  summation  by  parts 

on  the  RHS  of  (3.22)  yields 


E 


du+i  -  dv  d„  -  d„_ 


Ax, 


"  1  )  Ax 


|  nuhf 


Ax  ,  i 

^  ^  2 


E 

V 

y 

J  Ax 


(vu+ 1  -  v„,  d„+i  -  dv 


Av,,  ,  i ,  Ad, 


u+i 


E 


AA+i 


Av 


Ad^l  A 

■Av, 


2  ’  Av  1  UJr  2 


By  (3.23)  and  (3.24b),  the  semi-discrete  entropy  balance  amounts  to 
d 


dt 


y  u„(t)  aXv  =  ~y  — —  ( Av 


A 

2  Av, 


2  ’  Av  1  2 


(3.24a) 
.  (3.24b) 


(3.25) 


Here 


Ad, 


Av 


v+h  JZ=-\ 


dv(V,+I 


(0)  dt 


where  vv,  i(^)  is  given  by  (3.12).  By  the  admissibility  of  the  dissipative  NS  fluxes  dv  >  0  and  the 
RHS  of  (3.25)  is  indeed  non-positive.  Thus,  the  senri-discrete  scheme  (3.19a)  guarantees  the  total 
entropy  dissipation. 

In  the  specific  case  of  the  entropy  pair  (U,  F )  =  (— pS ,  — mS ),  the  entropy  variable  are  found  in 
(2.14),  and  we  explicitly  compute  the  inner  products  in  (3.24a)  as  (omitting  all  subscripts), 


E^Av.Ad) 


V 


The  discrete  entropy  balance  (1.11)  now  follows.  □ 

We  emphasize  the  main  point  made  here,  namely,  we  introduce  no  excessive  entropy  dissipation 
due  to  spurious,  artificial  numerical  viscosity:  by  (3.20),  the  senri-discrete  scheme  contains  the  precise 
amount  of  numerical  viscosity  to  enforce  the  correct  entropy  dissipation  dictated  by  the  NS  equations. 


3.3.  Time  discretization.  To  complete  the  computation  of  a  senri-discrete  scheme,  it  needs  to  be 
augmented  with  a  proper  time  discretization.  To  enable  a  large  tinre-stability  region  and  maintain 
simplicity,  the  three-stage  third-order  Runge-Kutta  (RK3)  method  will  be  used,  consult  [GST2001], 

uv(tn  +  At)  =  u„(t'1)  -| — —(211^1  +  3Kv>2  +  4 RE, 3), 

y 


(3.26) 
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where 


I<u,i 

Kv,  2 

K„,  3  = 


1 


Axy 

1 

Ax„ 


Ax, 


{r+,  (<,<+1) 


-f 


u 


u-li 


<) } 


n  At  At 

W  +  ~2~Kv,i,  u;/+i  + 

-C_1  (<_1  +  yif.-y.u"  +  ! 

3At  „  „  3At 


^+5 


_ ^ 


u;  + 


-Ku,2,  uv+1  + 


-K, 


v+1,2 


4  ■  4 

„  ,  3At  _  3At 

W-i  H  — H  -  Kv2 


The  resulting  fully-discrete  schemes  has  a  spatial  stencil  involving  seven-point  gridvalues,  with  two 
“ghost”  boundary  values  on  the  left  boundary  and  two  “ghost”  boundary  on  the  right  required  to  close 
the  system.  For  simplicity,  these  “ghost”  values  are  extrapolated  from  the  given  Dirichlet  boundary 
values.  We  note  in  passing  that  though  the  fully  explicit  RK3  time  discretization  need  not  conserve 
the  entropy,  it  introduces  a  negligible  amount  of  entropy  dissipation;  for  a  general  framework  of 
entropy-conservative  fully  discrete  schemes  consult  [LMR2002]. 


4.  Numerical  experiments 


We  consider  ideal  polytropic  gas  equations  as  an  approximation  of  air  with 
7  =  1.4,  Cv  =  716,  k  =  0.03,  X  +  2p  =  2.28  x  10“5 


We  simulate  the  Sod’s  shocktube  problem,  [Sodl978],  where  the  Euler  and  NS  equations  are  solved 
over  the  interval  [0,  1]  subject  to  Rienrann  initial  conditions 


(p,m,E)t= o  = 


(1.0,  0.0,  2.5)  0  <  x  <  0.5 

(0.125,  0.0,  0.25)  0.5  <  x  <  1. 


In  the  following  figures,  we  display  the  numerical  solutions  for  the  fully  discrete  scheme  (3.26)  with 
the  numerical  flux  (3.7),  or  in  its  equivalent  yet  simpler  form  (3.11).  Uniform  space  and  time  grid 
sizes,  Ax  and  At,  are  used.  Both  viscous  and  inviscid  cases  are  explored.  We  use  different  spatial 
resolutions  for  the  same  problem,  and  adjust  time  step  according  to  the  CFL  condition.  Different 
choices  of  entropy  function  are  also  tested  in  the  numerical  experiments.  We  group  our  results  into 
four  sets. 


1.  Euler  equations.  The  first  four  sets  of  figures  are  devoted  to  the  Euler  equations  with  zero 
viscous  and  heat  fluxes  (1.5). 

With  the  choice  of  the  entropy  pair 

(U( u),  F( u))  =  PP )T+^,  ^r^qipp)^^  ,  (4.1) 

V1  -  7  1  -  7  ) 

Figure  4.1  depicts  the  density,  velocity,  and  pressure  fields  at  t  =  0.05  and  t  =  0.1;  here  we  use 
Ax  =  0.001  and  ^  =  0.025.  Comparing  these  to  the  corresponding  results  of  the  canonical  entropy 
pair 

(U( u),  F( u))  =  (-pln(pp-'y),  — mln(p/T7))  ,  (4.2) 

in  figure  4.2,  we  see  that  the  different  choices  of  entropies  do  not  affect  the  behavior  of  the  numerical 
solutions.  Figures  4.1(d)  and  4.2(d)  demonstrate  the  conservation  of  the  total  entropies:  the  negligible 
amount  of  entropy  decay  ~  10~4  is  introduced  by  the  RK3  time  discretization. 

Next,  we  make  the  same  comparison  for  the  refined  the  spatial  mesh,  taking  Ax  =  0.00025,  ^  = 
0.1.  Figure  4.3  presents  the  computed  solutions  of  density,  velocity,  and  pressure  fields  at  t  =  0.05  and 
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t  =  0.1  with  the  entropy  pair  (4.1)  while  figure  4.4  depicts  the  solutions  with  the  canonical  entropy 
pair  (4.2).  The  total  entropy  is  shown  in  figures  4.3(d)  and  4.4(d). 

The  above  results  demonstrate  the  purely  dispersive  character  of  the  entropy  conservative  schemes. 
Dispersive  oscillations  on  the  mesh  scale  are  observed  in  shocks  and  contact  regions,  due  to  the  absence 
of  any  dissipation  mechanism,  consult  [Laxl986,  LevernroreLiul996].  The  numerical  solutions  do  not 
blow  up.  Actually,  as  we  refine  the  mesh,  these  dispersive  oscillations  approach  a  modulated  wave 
envelope.  The  study  of  these  modulated  waves  in  the  conservative  Euler  equations  would  be  an 
extremely  challenging  task.  A  similar  entropy  conservative  Lagrangian  formulation  of  Euler  equations 
of  [TrulioTriggerl961]  motivated  the  discussion  in  [Laxl986]. 

2  Navier-Stokes  equations  with  heat  flux.  We  solve  the  Navier-Stokes  equations  (1.4).  The 
results  are  summarized  in  the  next  three  sets  of  figures  4. 5-4. 7.  We  follow  the  same  pattern  of 
plotting  density,  velocity,  pressure  and  total  entropy.  As  before,  the  choice  of  entropy  pairs  (4.1)  in 
figures  4.5  and  4.6  are  very  similar. 

The  presence  of  heat  flux  causes  the  oscillations  to  be  dramatically  reduced  around  the  contact 
discontinuity.  Furthermore,  oscillations  are  significantly  damped  around  the  shock;  when  the  mesh 
is  well-refined,  figure  4.7  shows  that  heat  conduction  causes  these  oscillations  to  be  well  localized  in 
the  immediate  neighborhood  of  the  shocks.  If  the  mesh  is  underresolved,  a  small  portion  of  dispersive 
oscillations  persist  in  the  neighborhood  of  shocks. 

3.  Navier-stokes  equations  with  viscosity  and  no  heat  flux.  We  solve  the  viscous  NS  equations 
(1.3).  The  results  are  summarized  in  figures  4. 8-4. 9.  Since  the  results  are  essentially  independent  of 
the  choice  of  entropy,  we  chose  to  quote  here  only  the  results  for  the  canonical  pair  (4.2). 

The  viscosity  in  NS  equations  is  doing  a  better  job  than  heat  flux  in  damping  oscillations  around  the 
shock  discontinuity.  The  plots  of  total  entropy,  reveal  a  greater  entropy  decay  than  the  NS  equations 
with  heat  conduction.  On  the  other  hand,  we  still  observe  an  oscillatory  behavior  around  the  contact 
discontinuity,  even  with  the  refined  mesh  in  figure  4.9. 

4.  Full  Navier-stokes  equations  with  viscous  and  heat  fluxes.  In  figures  4.10-4.11  we  record 
the  results  for  the  full  NS  equations  (1.1).  As  before,  the  difference  due  to  different  entropy  functions 
is  undetectable  and  we  chose  to  record  here  only  the  canonical  entropy. 

As  expected,  these  numerical  solutions  are  the  smoothest  ones  found  in  our  numerical  experiments, 
especially  in  very  fine  meshes,  depicted  in  figure  4.11.  Small  oscillations  remain  with  underresolved 
meshes. 

Not  only  the  oscillations  around  the  shocks  are  damped  out  by  viscosity,  but  the  oscillations  around 
the  contact  discontinuity  are  significantly  reduced  due  to  the  heat  flux.  Compared  with  the  results 
of  NS  equations  with  heat  conduction  (1.4)  in  figures  4. 6-4. 7,  oscillations  in  the  neighborhood  of  the 
shock  are  better  damped  here  thanks  to  the  viscosity  terms.  The  remaining  sharp  “spike”  at  the  tip 
of  shock  discontinuity  is  due  to  the  relatively  small  viscosity  coefficient  of  air. 
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Figure  4.1.  Euler  equations  with  1000  spatial  gridpoints,  U( u)  =  -(pp) 1+7  ,  A t  = 

2.5  x  10“5,  Ax  =  10“3 
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1  -D  Euler, velocity ji=1 ,4,CV=71 6,K=0,entropy=-p  ln(pp-1),At=2.5e-005,Ax=0.001  ,Tm  =0.1 


(a)  Density 


(b)  Velocity 


1-D  Euler, pressure, y=1. 4, Cv=71 6, K=0,entropy=-p  ln(pp  'l),At=2.5e-005,Ax=0.001,T  =0.1 


(c)  Pressure 


(d)  Total  entropy  v.s.  time 


Figure  4.2.  Euler  equations  with  1000  spatial  gridpoints,  U( u)  =  —pln(pp  7)  and 
same  At  and  Ax  as  figure  4.1 
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1-D  Euler,density,Y=1.4,Cv=716,K=0,entropy=(1+^/(1-^‘(PP)1/(1^=2-5e-005,Ax=0.00025,Tmax=0.1  i-D  Euler,velocityjfc1.4,Cv=716,K=0,entropy=(1+^/(1-^*(p  p)1/<1+^At=2.5e-005,Ax=0.00025,Tmax=0. 
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1  -D  Euler, pressure, y=1 .4,Cy-71 6,K-0,entropy-(1  +^/(1  -y)*(p  p)1/^1+^,At-2.5e-005, Ax-0.00025, Tmax-0  i  _q  Euler, entropy-conservative  scheme  w/  3rd  R-K  mtd, entropy  vs  time, entropy  =  (1  -fy)/(1  -^*(p  p)1^ 


Figure  4.3.  Euler  equations  with  4000  spatial  gridpoints,  U{ u)  =  -{pp) 1+7  ,  At  = 

2.5  x  10-5,  Ax  =  2.5  x  10"4 
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1-D  Euler, density, y=1. 4, C  =71 6,K=0,entropy=-p  ln(pp  7),At=2.5e-005,Ax=0.00025,T  =0.1 


1-D  Euler,velocityy=1.4,Cv=716,K=0,entropy=-p  ln(pp  ^),At=2.5e-005,Ax=0.00025,Tmax=0.1 
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(d)  Total  entropy  v.s.  time 


Figure  4.4.  Euler  equations  with  4000  spatial  gridpoints,  U( u)  =  —  pln(pp  7)  and 
same  At  and  Ax  as  Figure  4.3 
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(a)  Density 


(b)  Velocity 


Figure  4.5.  Navier-Stokes  equations  with  heat  conduction  and  no  viscous  term.  1000 
spatial  gridpoints,  U(u)  =  ■  (pp)~ ,  A t  =  2.5  x  10~5,  Ax  =  10-3 
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1-D  Euler, density;^  1 .4,Cv=716,K=0.03,entropy=-p  ln(pp  i),At=2.5e-005,Ax=0.001  -T  =0.1  1-D  Euler, velocity;^  1 .4,Cv=716,K=0.03,entropy=-p  ln(pp  ^,At=2.5e-005,Ax=0.001  ,Tm  =0.1 


(a)  Density 


(b)  Velocity 


1-D  Euler, pressure, y=1 .4, Cv=71 6, K=0.03,entropy=-p  ln(pp  Y),At=2.5e-005,Ax=0.001,Tmax=0.1  i-D  Euler,  entropy-conservative  scheme  w /  3rd  R-K  mtd,  entropy  vs  time, entropy  =  -p  ln(pp  1) 


(c)  Pressure 


(d)  Total  entropy  v.s.  time 


Figure  4.6.  Navier-Stokes  equations  with  heat  conduction  and  no  viscous  terms.  1000 
spatial  gridpoints,  U( u)  =  — pln(p/9~7)  and  same  At  and  Ax  as  Figure  4.5 
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1-D  Euler,densityy=1.4,Cv=716,K=0.03,entropy=-p  ln(pp  1),At=2.5e-005,Ax=0.00025,Tmax=0.1  1-D  Euler, velocity;^ ,4,CV=71 6,K=0.03,entropy=-p  ln(pp-1),At=2.5e-005,Ax=0.00025,Tmax=0.1 


(a)  Density 


(b)  Velocity 


(c)  Pressure 


(d)  Total  entropy  v.s.  time 


Figure  4.7.  Navier-Stokes  equations  with  heat  conduction  and  no  viscous  terms.  4000 
spatial  grids,  U( u)  =  —  pln(pp~~l) ,  At  =  2.5  x  10~5,  Ax  =  2.5  x  10-4 
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1-D  N-S,density1A.+2p=2.28e-005,y=1.4,Cv=716,K=0,entropy=-p  ln(pp  ^),At=2.5e-005,Ax=0.001  1-D  N-S,  velocity  A+2p=2.28e-005,y=1. 4, Cv=71 6, K=0,entropy=-p  ln(pp  'l),At=2.5e-005,Ax=0.001 


(a)  Density 


(b)  Velocity 


1  -D  N-S, pressure, A.+2p=2.28e-005,,y=1 .4,CV=71 6,K=0,entropy=-p  ln(pp-^),At=2.5e-005,Ax=0.001  i  -D  N-S,  entropy-conservative  scheme  w/  3rd  R-K  mtd, entropy  vs  time, entropy  =  -p  ln(pp  ^ 


(c)  Pressure 


(d)  Total  entropy  v.s.  time 


Figure  4.8.  Navier-Stokes  equations  with  viscous  terms  and  no  heat  conduction.  1000 
spatial  gridpoints,  U( u)  =  —  pln(pp~1) ,  At  =  2.5  x  10~5,  Ax  =  10~3 
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1  -D  N-S, density A+2(i=2.28e-005, 7^1 .4,CV=71 6,K=0,entropy=-p  ln(pp_Y),At=2.5e-005,Ax=0.00025  1  -D  N-S, velocity  A+2p=2.28e-005,7fc1 ,4,CV=71 6,K=0,entropy=-p  ln(pp  7),At=2.5e-005,Ax=0.00025 


(a)  Density 


(b)  Velocity 


Figure  4.9.  Navier-Stokes  equations  with  viscous  terms  and  no  heat  conduction.  4000 
spatial  grids,  U( u)  =  —  pln(pp~'y) ,  At  =  2.5  x  10~5,  Ax  =  2.5  x  10-4 
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1-D  N-S,density,A.+2p=2.28e-005,Y=1.4,Cv=716,K=0.03,entropy=-p  ln(pp  ^,At=2.5e-005,Ax=0.001  1-D  N-S, velocity A+2p=2.28e-005,,y=1. 4, Cv=71 6, K=0.03,entropy=-p  ln(pp  ^),At=2.5e-005,Ax=0.001 


(a)  Density 


(b)  Velocity 


1  -D  N-S, pressure, A.+2p=2.28e-005,Y=1 .4,C  =71 6,K=0.03,entropy=-p  ln(pp  ^),At=2.5e-005,Ax=0.00  i  -D  N-S,  entropy-conservative  scheme  w/  3rd  R-K  mtd, entropy  vs  time, entropy  =  -p  ln(pp-1) 


(c)  Pressure 


(d)  Total  entropy  v.s.  time 


Figure  4.10.  Navier-Stokes  equations  with  viscosity  and  heat  conduction.  1000  spa¬ 
tial  gridpoints,  U{ u)  =  —  pln(pp~'y) ,  At  =  2.5  x  10”5,  Ax  =  10~3 
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1-D  N-S,densityA+2|i=2.28e-005,7!=1.4,Cv=716,K=0.03,entropy=-p  ln(pp  ^),At=2.5e-005,Ax=0.0002  1-D  N-S,velocityA+2p=2.28e-005,y=1.4,Cv=716,K=0.03,entropy=-p  ln(pp  'l),At=2.5e-005,Ax=0.0002 


(a)  Density 


(b)  Velocity 


1-D  N-S, pressure, A.+2p=2.28e-005,y=1. 4, Cy=71 6, K=0.03,entropy=-p  ln(pp  1),At=2.5e-005,Ax=0.000; 


(c)  Pressure 


(d)  Total  entropy  v.s.  time 


Figure  4.11.  Navier-Stokes  equations  with  viscosity  and  heat  conduction.  4000  spa¬ 
tial  gridpoints,  U( u)  =  —pln(pp~'y) ,  At  =  2.5  x  10”5,  Ax  =  2.5  x  10~4 
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